knitr::opts_chunk$set(echo = TRUE)
knitr::opts_knit$set(root.dir = "~/DADS")
library(phyloseq); packageVersion("phyloseq")
## [1] '1.50.0'
library(ggplot2); packageVersion("ggplot2")
## [1] '3.5.1'
library(readr); packageVersion("readr")
## [1] '2.1.5'
library(purrr); packageVersion("purrr")
## [1] '1.0.4'
library(furrr); packageVersion("furrr")
## Loading required package: future
## [1] '0.3.1'
library(dplyr); packageVersion("dplyr")
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
## [1] '1.1.4'
library(stringr); packageVersion("stringr")
## [1] '1.5.1'
library(metacoder); packageVersion("metacoder")
## This is metacoder version 0.3.8 (stable)
##
## Attaching package: 'metacoder'
## The following object is masked from 'package:ggplot2':
##
## map_data
## The following object is masked from 'package:phyloseq':
##
## filter_taxa
## [1] '0.3.8'
library(data.table); packageVersion("data.table")
##
## Attaching package: 'data.table'
## The following objects are masked from 'package:dplyr':
##
## between, first, last
## The following object is masked from 'package:purrr':
##
## transpose
## [1] '1.16.4'
library(Biostrings); packageVersion("Biostrings")
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following objects are masked from 'package:dplyr':
##
## combine, intersect, setdiff, union
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, aperm, append, as.data.frame, basename, cbind,
## colnames, dirname, do.call, duplicated, eval, evalq, Filter, Find,
## get, grep, grepl, intersect, is.unsorted, lapply, Map, mapply,
## match, mget, order, paste, pmax, pmax.int, pmin, pmin.int,
## Position, rank, rbind, Reduce, rownames, sapply, saveRDS, setdiff,
## table, tapply, union, unique, unsplit, which.max, which.min
## Loading required package: S4Vectors
## Loading required package: stats4
##
## Attaching package: 'S4Vectors'
## The following objects are masked from 'package:data.table':
##
## first, second
## The following objects are masked from 'package:dplyr':
##
## first, rename
## The following object is masked from 'package:utils':
##
## findMatches
## The following objects are masked from 'package:base':
##
## expand.grid, I, unname
## Loading required package: IRanges
##
## Attaching package: 'IRanges'
## The following object is masked from 'package:data.table':
##
## shift
## The following object is masked from 'package:metacoder':
##
## reverse
## The following objects are masked from 'package:dplyr':
##
## collapse, desc, slice
## The following object is masked from 'package:purrr':
##
## reduce
## The following object is masked from 'package:phyloseq':
##
## distance
## Loading required package: XVector
##
## Attaching package: 'XVector'
## The following object is masked from 'package:purrr':
##
## compact
## Loading required package: GenomeInfoDb
##
## Attaching package: 'Biostrings'
## The following object is masked from 'package:metacoder':
##
## complement
## The following object is masked from 'package:base':
##
## strsplit
## [1] '2.74.1'
library(vegan);packageVersion("vegan")
## Loading required package: permute
## Loading required package: lattice
## [1] '2.6.10'
library(agricolae);packageVersion("agricolae")
## [1] '1.3.7'
library(magrittr);packageVersion("magrittr")
##
## Attaching package: 'magrittr'
## The following object is masked from 'package:purrr':
##
## set_names
## [1] '2.0.3'
library(data.table);packageVersion("data.table")
## [1] '1.16.4'
library(DESeq2);packageVersion("DESeq2")
## Loading required package: GenomicRanges
##
## Attaching package: 'GenomicRanges'
## The following object is masked from 'package:magrittr':
##
## subtract
## Loading required package: SummarizedExperiment
## Loading required package: MatrixGenerics
## Loading required package: matrixStats
##
## Attaching package: 'matrixStats'
## The following object is masked from 'package:dplyr':
##
## count
##
## Attaching package: 'MatrixGenerics'
## The following objects are masked from 'package:matrixStats':
##
## colAlls, colAnyNAs, colAnys, colAvgsPerRowSet, colCollapse,
## colCounts, colCummaxs, colCummins, colCumprods, colCumsums,
## colDiffs, colIQRDiffs, colIQRs, colLogSumExps, colMadDiffs,
## colMads, colMaxs, colMeans2, colMedians, colMins, colOrderStats,
## colProds, colQuantiles, colRanges, colRanks, colSdDiffs, colSds,
## colSums2, colTabulates, colVarDiffs, colVars, colWeightedMads,
## colWeightedMeans, colWeightedMedians, colWeightedSds,
## colWeightedVars, rowAlls, rowAnyNAs, rowAnys, rowAvgsPerColSet,
## rowCollapse, rowCounts, rowCummaxs, rowCummins, rowCumprods,
## rowCumsums, rowDiffs, rowIQRDiffs, rowIQRs, rowLogSumExps,
## rowMadDiffs, rowMads, rowMaxs, rowMeans2, rowMedians, rowMins,
## rowOrderStats, rowProds, rowQuantiles, rowRanges, rowRanks,
## rowSdDiffs, rowSds, rowSums2, rowTabulates, rowVarDiffs, rowVars,
## rowWeightedMads, rowWeightedMeans, rowWeightedMedians,
## rowWeightedSds, rowWeightedVars
## Loading required package: Biobase
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Attaching package: 'Biobase'
## The following object is masked from 'package:MatrixGenerics':
##
## rowMedians
## The following objects are masked from 'package:matrixStats':
##
## anyMissing, rowMedians
## The following object is masked from 'package:phyloseq':
##
## sampleNames
## [1] '1.46.0'
library(cowplot);packageVersion("cowplot")
## [1] '1.1.3'
library(doParallel);packageVersion("doParallel")
## Loading required package: foreach
##
## Attaching package: 'foreach'
## The following objects are masked from 'package:purrr':
##
## accumulate, when
## Loading required package: iterators
## Loading required package: parallel
## [1] '1.0.17'
Load data from all three experiments
load("Data/ITS/filtered_data/dadsexp.RData")
dads_its.ps<-dads.ps
load("Data/ITS/filtered_data/dads_5asv.RData")
dads_5asv_its.ps<-dads_5asv.ps
load("Data/16S/filtered_data/dadsexp.RData")
dads_16S.ps<-dads.ps
load("Data/16S/filtered_data/dads_5asv.RData")
dads_5asv_16S.ps<-dads_5asv.ps
load("Data/RPS10/filtered_data/dadsexp.RData")
dads_rps10.ps<-dads.ps
load("Data/RPS10/filtered_data/dads_5asv.RData")
dads_5asv_rps10.ps<-dads_5asv.ps
Let’s calculate alpha diversity metrics for ITS dataset
set.seed(1)
variables <- c("Treatment", "Soil", "Incubation")
measures <- c("Observed", "InvSimpson")
AD_its <- phyloseq::estimate_richness(dads_its.ps, measures = c("Observed", "InvSimpson"))
AD_its <- cbind(AD_its, sample_data(dads_its.ps)[rownames(AD_its), c("Treatment", "Soil", "Incubation")])
variables <- c("Treatment", "Soil", "Incubation")
plots_its <- list()
for (variable in variables) {
p <- phyloseq::plot_richness(dads_its.ps, x = variable, measures = c("Observed", "InvSimpson")) +
geom_boxplot() +
theme_classic() + # Use theme_classic() as a starting point
theme(
panel.border = element_rect(color = "black", fill = NA, size = 1), # Add border around plot panel
axis.line = element_blank(), # Remove axis lines (they'll be replaced by the panel border)
plot.margin = unit(c(1, 1, 1, 1), "lines"), # Add some margin around the plot
plot.title = element_text(margin = margin(b = 10)), # Add some space below the title
axis.title.x = element_text(margin = margin(t = 10)), # Add some space above the x-axis title
axis.title.y = element_text(margin = margin(r = 10)) # Add some space to the right of the y-axis title
)
plots_its[[variable]] <- p
}
## Warning: The `size` argument of `element_rect()` is deprecated as of ggplot2 3.4.0.
## ℹ Please use the `linewidth` argument instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
#Let's calculate Tukey results
# Function to perform Tukey HSD and extract letters
get_tukey_letters <- function(data, measure, variable) {
formula <- as.formula(paste(measure, "~", variable))
aov_result <- aov(formula, data = data)
tukey_result <- HSD.test(aov_result, variable, group = TRUE)
letters <- tukey_result$groups$groups
return(letters)
}
# Perform Tukey HSD tests and store results
tukey_results_ITS <- list()
for (variable in variables) {
for (measure in measures) {
tryCatch({
letters <- get_tukey_letters(AD_its, measure, variable)
tukey_results_ITS[[paste(measure, variable, sep = "_")]] <- letters
}, error = function(e) {
message("Error in Tukey test for ", measure, " by ", variable, ": ", e$message)
})
}
}
print(tukey_results_ITS)
## $Observed_Treatment
## [1] "a" "b"
##
## $InvSimpson_Treatment
## [1] "a" "b"
##
## $Observed_Soil
## [1] "a" "a"
##
## $InvSimpson_Soil
## [1] "a" "a"
##
## $Observed_Incubation
## [1] "a" "a"
##
## $InvSimpson_Incubation
## [1] "a" "a"
Let’s calculate alpha diversity metrics for 16S dataset
set.seed(1)
variables <- c("Treatment", "Soil", "Incubation")
measures <- c("Observed", "InvSimpson")
AD_16S <- phyloseq::estimate_richness(dads_16S.ps, measures = c("Observed", "InvSimpson"))
rownames(AD_16S) <- gsub("\\X", "", rownames(AD_16S))
# Combine richness estimates with sample metadata
AD_16S <- cbind(AD_16S, sample_data(dads_16S.ps)[rownames(AD_16S), c("Treatment", "Soil", "Incubation")])
variables <- c("Treatment", "Soil", "Incubation")
# Create an empty list to store plots
plots_16S <- list()
# Loop through each variable and create the plot
for (variable in variables) {
p <- phyloseq::plot_richness(dads_16S.ps, x = variable, measures = c("Observed", "InvSimpson")) +
geom_boxplot() +
theme_classic() + # Use theme_classic() as a starting point
theme(
panel.border = element_rect(color = "black", fill = NA, size = 1), # Add border around plot panel
axis.line = element_blank(), # Remove axis lines (they'll be replaced by the panel border)
plot.margin = unit(c(1, 1, 1, 1), "lines"), # Add some margin around the plot
plot.title = element_text(margin = margin(b = 10)), # Add some space below the title
axis.title.x = element_text(margin = margin(t = 10)), # Add some space above the x-axis title
axis.title.y = element_text(margin = margin(r = 10)) # Add some space to the right of the y-axis title
)
plots_16S[[variable]] <- p
}
#Let's calculate Tukey results
# Function to perform Tukey HSD and extract letters
get_tukey_letters <- function(data, measure, variable) {
formula <- as.formula(paste(measure, "~", variable))
aov_result <- aov(formula, data = data)
tukey_result <- HSD.test(aov_result, variable, group = TRUE)
letters <- tukey_result$groups$groups
return(letters)
}
# Perform Tukey HSD tests and store results
tukey_results_16S <- list()
for (variable in variables) {
for (measure in measures) {
tryCatch({
letters <- get_tukey_letters(AD_16S, measure, variable)
tukey_results_16S[[paste(measure, variable, sep = "_")]] <- letters
}, error = function(e) {
message("Error in Tukey test for ", measure, " by ", variable, ": ", e$message)
})
}
}
print(tukey_results_16S)
## $Observed_Treatment
## [1] "a" "b"
##
## $InvSimpson_Treatment
## [1] "a" "b"
##
## $Observed_Soil
## [1] "a" "a"
##
## $InvSimpson_Soil
## [1] "a" "a"
##
## $Observed_Incubation
## [1] "a" "a"
##
## $InvSimpson_Incubation
## [1] "a" "a"
Let’s calculate alpha diversity metrics for RPS10 dataset
set.seed(1)
variables <- c("Treatment", "Soil", "Incubation")
measures <- c("Observed", "InvSimpson")
# Estimate alpha diversity
AD_RPS10 <- estimate_richness(dads_rps10.ps, measures = c("Observed", "InvSimpson"))
## Warning in estimate_richness(dads_rps10.ps, measures = c("Observed", "InvSimpson")): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.
# Combine richness estimates with sample metadata
AD_RPS10 <- cbind(AD_RPS10, sample_data(dads_rps10.ps)[rownames(AD_RPS10), c("Treatment", "Soil", "Incubation")])
variables <- c("Treatment", "Soil", "Incubation")
# Create an empty list to store plots
plots_rps10 <- list()
# Loop through each variable and create the plot
for (variable in variables) {
p <- phyloseq::plot_richness(dads_rps10.ps, x = variable, measures = c("Observed", "InvSimpson")) +
geom_boxplot() +
theme_classic() + # Use theme_classic() as a starting point
theme(
panel.border = element_rect(color = "black", fill = NA, size = 1), # Add border around plot panel
axis.line = element_blank(), # Remove axis lines (they'll be replaced by the panel border)
plot.margin = unit(c(1, 1, 1, 1), "lines"), # Add some margin around the plot
plot.title = element_text(margin = margin(b = 10)), # Add some space below the title
axis.title.x = element_text(margin = margin(t = 10)), # Add some space above the x-axis title
axis.title.y = element_text(margin = margin(r = 10)) # Add some space to the right of the y-axis title
)
plots_rps10[[variable]] <- p
}
## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.
## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.
## Warning in estimate_richness(physeq, split = TRUE, measures = measures): The data you have provided does not have
## any singletons. This is highly suspicious. Results of richness
## estimates (for example) are probably unreliable, or wrong, if you have already
## trimmed low-abundance taxa from the data.
##
## We recommended that you find the un-trimmed data and retry.
#Let's calculate Tukey results
# Function to perform Tukey HSD and extract letters
get_tukey_letters <- function(data, measure, variable) {
formula <- as.formula(paste(measure, "~", variable))
aov_result <- aov(formula, data = data)
tukey_result <- HSD.test(aov_result, variable, group = TRUE)
letters <- tukey_result$groups$groups
return(letters)
}
# Perform Tukey HSD tests and store results
tukey_results_RPS10 <- list()
for (variable in variables) {
for (measure in measures) {
tryCatch({
letters <- get_tukey_letters(AD_RPS10, measure, variable)
tukey_results_RPS10[[paste(measure, variable, sep = "_")]] <- letters
}, error = function(e) {
message("Error in Tukey test for ", measure, " by ", variable, ": ", e$message)
})
}
}
print(tukey_results_RPS10)
## $Observed_Treatment
## [1] "a" "a"
##
## $InvSimpson_Treatment
## [1] "a" "a"
##
## $Observed_Soil
## [1] "a" "a"
##
## $InvSimpson_Soil
## [1] "a" "a"
##
## $Observed_Incubation
## [1] "a" "a"
##
## $InvSimpson_Incubation
## [1] "a" "a"
Let’s combine the treatment plots
create_label <- function(text, x = 0.5, y = 0.95) { # Centered the label horizontally
ggdraw() +
draw_label(text, fontface = 'bold', x = x, y = y, hjust = 0.5, vjust = 1) # Changed hjust to 0.5 for center alignment
}
label_A <- create_label("A.")
label_B <- create_label("B.")
label_C <- create_label("C.")
# Arrange plots
plots_arrange_its <- gridExtra::grid.arrange(
plots_its[["Treatment"]]
)
plots_arrange_16S <- gridExtra::grid.arrange(
plots_16S[["Treatment"]]
)
plots_arrange_rps10 <- gridExtra::grid.arrange(
plots_rps10[["Treatment"]]
)
combined_plots_16S <- plot_grid(label_A, plots_arrange_16S, ncol = 2, rel_widths = c(0.1, 1)) # Reduced label width
combined_plots_its <- plot_grid(label_B, plots_arrange_its, ncol = 2, rel_widths = c(0.1, 1)) # Reduced label width
combined_plots_rps10 <- plot_grid(label_C, plots_arrange_rps10, ncol = 2, rel_widths = c(0.1, 1)) # Reduced label width
facetted_figure_trt <- plot_grid(
combined_plots_16S,
combined_plots_its,
combined_plots_rps10,
nrow = 3,
align = 'v',
axis = 'l'
)
facetted_figure_trt_with_margins <- plot_grid(
NULL, facetted_figure_trt, NULL,
ncol = 3,
rel_widths = c(0, 1, 0.08) # 5% margin on each side
)
facetted_figure_trt_with_margins
# Save the figure
ggsave("Figures/facetted_figure_alphadiv_treatment.svg", facetted_figure_trt_with_margins, width=8, height=12)
ggsave("Figures/facetted_figure_alphadiv_treatment.pdf", facetted_figure_trt_with_margins, width=8, height=12)
Let’s combine the soil type plots
create_label <- function(text, x = 0.5, y = 0.95) { # Centered the label horizontally
ggdraw() +
draw_label(text, fontface = 'bold', x = x, y = y, hjust = 0.5, vjust = 1) # Changed hjust to 0.5 for center alignment
}
label_A <- create_label("A.")
label_B <- create_label("B.")
label_C <- create_label("C.")
# Arrange plots
plots_arrange_its <- gridExtra::grid.arrange(
plots_its[["Soil"]]
)
plots_arrange_16S <- gridExtra::grid.arrange(
plots_16S[["Soil"]]
)
plots_arrange_rps10 <- gridExtra::grid.arrange(
plots_rps10[["Soil"]]
)
combined_plots_16S <- plot_grid(label_A, plots_arrange_16S, ncol = 2, rel_widths = c(0.1, 1)) # Reduced label width
combined_plots_its <- plot_grid(label_B, plots_arrange_its, ncol = 2, rel_widths = c(0.1, 1)) # Reduced label width
combined_plots_rps10 <- plot_grid(label_C, plots_arrange_rps10, ncol = 2, rel_widths = c(0.1, 1)) # Reduced label width
facetted_figure_soil <- plot_grid(
combined_plots_16S,
combined_plots_its,
combined_plots_rps10,
nrow = 3,
align = 'v',
axis = 'l'
)
facetted_figure_soil_with_margins <- plot_grid(
NULL, facetted_figure_soil, NULL,
ncol = 3,
rel_widths = c(0, 1, 0.08) # 5% margin on each side
)
facetted_figure_soil_with_margins
# Save the figure
ggsave("Figures/facetted_figure_alphadiv_soil.svg", facetted_figure_soil_with_margins, width=8, height=12)
ggsave("Figures/facetted_figure_alphadiv_soil.pdf", facetted_figure_soil_with_margins, width=8, height=12)
Let’s combine the incubation plots
create_label <- function(text, x = 0.5, y = 0.95) { # Centered the label horizontally
ggdraw() +
draw_label(text, fontface = 'bold', x = x, y = y, hjust = 0.5, vjust = 1) # Changed hjust to 0.5 for center alignment
}
label_A <- create_label("A.")
label_B <- create_label("B.")
label_C <- create_label("C.")
# Arrange plots
plots_arrange_its <- gridExtra::grid.arrange(
plots_its[["Incubation"]]
)
plots_arrange_16S <- gridExtra::grid.arrange(
plots_16S[["Incubation"]]
)
plots_arrange_rps10 <- gridExtra::grid.arrange(
plots_rps10[["Incubation"]]
)
combined_plots_16S <- plot_grid(label_A, plots_arrange_16S, ncol = 2, rel_widths = c(0.1, 1)) # Reduced label width
combined_plots_its <- plot_grid(label_B, plots_arrange_its, ncol = 2, rel_widths = c(0.1, 1)) # Reduced label width
combined_plots_rps10 <- plot_grid(label_C, plots_arrange_rps10, ncol = 2, rel_widths = c(0.1, 1)) # Reduced label width
facetted_figure_inc <- plot_grid(
combined_plots_16S,
combined_plots_its,
combined_plots_rps10,
nrow = 3,
align = 'v',
axis = 'l'
)
facetted_figure_inc_with_margins <- plot_grid(
NULL, facetted_figure_inc, NULL,
ncol = 3,
rel_widths = c(0, 1, 0.08) # 5% margin on each side
)
facetted_figure_inc_with_margins
# Save the figure
ggsave("Figures/facetted_figure_alphadiv_incubation.svg", facetted_figure_inc_with_margins, width=8, height=12)
ggsave("Figures/facetted_figure_alphadiv_incubation.pdf", facetted_figure_inc_with_margins, width=8, height=12)
Let’s proceed with ordination for the ITS dataset-PCoA
set.seed(1)
strat_rel_dads_its.ps<-transform_sample_counts(dads_5asv_its.ps, function(x) x / sum(x))
ps_PCoA_its = ordinate(strat_rel_dads_its.ps,
method="PCoA",
distance="bray")
plot_its_trt<- plot_ordination(strat_rel_dads_its.ps,
ps_PCoA_its ,
color= "Treatment") +
#ggtitle("PCoA ordination plot - relative abundance") +
labs(color = "Treatment") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_its_trt<-plot_its_trt +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_its_trt
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
plot_its_inc<- plot_ordination(strat_rel_dads_its.ps,
ps_PCoA_its ,
color= "Incubation") +
#ggtitle("PCoA ordination plot - relative abundance") +
labs(color = "Incubation") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_its_inc<-plot_its_inc +
stat_ellipse(linetype =2) +
theme_bw()
plot_ellipses_its_inc
plot_its_soil<- plot_ordination(strat_rel_dads_its.ps,
ps_PCoA_its ,
color= "Soil") +
#ggtitle("PCoA ordination plot - relative abundance") +
labs(color = "Soil type") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_its_soil<-plot_its_soil +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_its_soil
Let’s proceed with ordination for the 16S dataset-PCoA
set.seed(1)
strat_rel_dads_16S.ps<-transform_sample_counts(dads_5asv_16S.ps, function(x) x / sum(x))
ps_PCoA_16S = ordinate(strat_rel_dads_16S.ps,
method="PCoA",
distance="bray")
plot_16S_trt<- plot_ordination(strat_rel_dads_16S.ps,
ps_PCoA_16S ,
color= "Treatment") +
#ggtitle("PCoA ordination plot - relative abundance") +
labs(color = "Treatment") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_16S_trt<-plot_16S_trt +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_16S_trt
plot_16S_inc<- plot_ordination(strat_rel_dads_16S.ps,
ps_PCoA_16S ,
color= "Incubation") +
#ggtitle("PCoA ordination plot - relative abundance") +
labs(color = "Incubation") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_16S_inc<-plot_16S_inc +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_16S_inc
plot_16S_soil<- plot_ordination(strat_rel_dads_16S.ps,
ps_PCoA_16S ,
color= "Soil") +
#ggtitle("PCoA ordination plot - relative abundance") +
labs(color = "Soil type") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_16S_soil<-plot_16S_soil +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_16S_soil
Let’s proceed with ordination for the RPS10 dataset-pCoA
set.seed(1)
strat_rel_dads_rps10.ps<-transform_sample_counts(dads_5asv_rps10.ps, function(x) x / sum(x))
ps_PCoA_rps10 = ordinate(strat_rel_dads_rps10.ps,
method="PCoA",
distance="bray")
plot_rps10_trt<- plot_ordination(strat_rel_dads_rps10.ps,
ps_PCoA_rps10 ,
color= "Treatment") +
#ggtitle("PCoA ordination plot - relative abundance") +
labs(color = "Treatment") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_rps10_trt<-plot_rps10_trt +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_rps10_trt
plot_rps10_inc<- plot_ordination(strat_rel_dads_rps10.ps,
ps_PCoA_rps10 ,
color= "Incubation") +
#ggtitle("PCoA ordination plot - relative abundance") +
labs(color = "Incubation") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_rps10_inc<-plot_rps10_inc +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_rps10_inc
plot_rps10_soil<- plot_ordination(strat_rel_dads_rps10.ps,
ps_PCoA_rps10 ,
color= "Soil") +
#ggtitle("PCoA ordination plot - relative abundance") +
labs(color = "Soil type") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_rps10_soil<-plot_rps10_soil +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_rps10_soil
Let’s combine plots into facetted plot
### Let's combine plots into facetted plot
create_column <- function(ord_plot) {
plot_grid(
ord_plot,
nrow = 1,
align = 'v',
axis = 'l',
rel_heights = c(1)
)
}
# Create columns for each dataset and factor
col_16S_treatment <- create_column(plot_ellipses_16S_trt)
col_16S_soil <- create_column(plot_ellipses_16S_soil)
col_16S_incubation <- create_column(plot_ellipses_16S_inc)
col_its_treatment <- create_column(plot_ellipses_its_trt)
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
col_its_soil <- create_column(plot_ellipses_its_soil)
col_its_incubation <- create_column(plot_ellipses_its_inc)
col_rps10_treatment <- create_column(plot_ellipses_rps10_trt)
col_rps10_soil <- create_column(plot_ellipses_rps10_soil)
col_rps10_incubation <- create_column(plot_ellipses_rps10_inc)
# Function to create headers with horizontal lines above
create_header <- function(text) {
ggplot() +
geom_hline(yintercept = 0.8, color = "black", size = 0.5) +
geom_text(aes(x = 0.02, y = 0.4, label = text), size = 5, fontface = "bold", hjust = 0) +
theme_void() +
theme(plot.margin = margin(5, 0, 0, 10, "pt")) +
xlim(0, 1) + ylim(0, 1)
}
# Create row plots with labels
row_16S <- plot_grid(col_16S_treatment, col_16S_soil, col_16S_incubation,
ncol = 3, labels = "A.", label_size = 14, hjust = -0.5, vjust = 2)
row_its <- plot_grid(col_its_treatment, col_its_soil, col_its_incubation,
ncol = 3, labels = "B.", label_size = 14, hjust = -0.5, vjust = 2)
row_rps10 <- plot_grid(col_rps10_treatment, col_rps10_soil, col_rps10_incubation,
ncol = 3, labels = "C.", label_size = 14, hjust = -0.5, vjust = 2)
# Combine all elements
facetted_figure <- plot_grid(
row_16S,
row_its,
row_rps10,
nrow = 3,
align = 'v',
axis = 'l',
rel_heights = c(1.1, 1, 1) # Slightly more height for the first row
)
facetted_figure
# Save the figure
ggsave("Figures/facetted_figure_with_ordination_PCoA.svg", facetted_figure, width = 12, height = 12)
ggsave("Figures/facetted_figure_with_ordination_PCoA.pdf", facetted_figure, width = 12, height = 12)
Let’s proceed with ordination for the ITS dataset-NMDS
set.seed(1)
strat_rel_dads_its.ps<-transform_sample_counts(dads_5asv_its.ps, function(x) x / sum(x))
ps_NMDS_its = ordinate(strat_rel_dads_its.ps,
method="NMDS",
distance="bray")
## Run 0 stress 9.209135e-05
## Run 1 stress 9.711561e-05
## ... Procrustes: rmse 4.557734e-05 max resid 0.000145828
## ... Similar to previous best
## Run 2 stress 9.941333e-05
## ... Procrustes: rmse 4.783054e-05 max resid 0.0001492673
## ... Similar to previous best
## Run 3 stress 9.765011e-05
## ... Procrustes: rmse 5.327067e-05 max resid 0.0001430395
## ... Similar to previous best
## Run 4 stress 9.298127e-05
## ... Procrustes: rmse 9.231115e-05 max resid 0.0002031309
## ... Similar to previous best
## Run 5 stress 9.956743e-05
## ... Procrustes: rmse 6.453389e-05 max resid 0.0001518031
## ... Similar to previous best
## Run 6 stress 9.81565e-05
## ... Procrustes: rmse 6.29341e-05 max resid 0.0001313132
## ... Similar to previous best
## Run 7 stress 9.415583e-05
## ... Procrustes: rmse 2.284447e-05 max resid 8.763404e-05
## ... Similar to previous best
## Run 8 stress 9.704892e-05
## ... Procrustes: rmse 0.0001012395 max resid 0.0002088867
## ... Similar to previous best
## Run 9 stress 9.824374e-05
## ... Procrustes: rmse 6.666882e-05 max resid 0.0001454821
## ... Similar to previous best
## Run 10 stress 9.985598e-05
## ... Procrustes: rmse 0.000100697 max resid 0.0002126261
## ... Similar to previous best
## Run 11 stress 9.996398e-05
## ... Procrustes: rmse 6.547809e-05 max resid 0.0001601838
## ... Similar to previous best
## Run 12 stress 9.761524e-05
## ... Procrustes: rmse 4.580355e-05 max resid 0.0001486625
## ... Similar to previous best
## Run 13 stress 9.698329e-05
## ... Procrustes: rmse 5.20759e-05 max resid 0.0001371574
## ... Similar to previous best
## Run 14 stress 9.946996e-05
## ... Procrustes: rmse 5.137318e-05 max resid 0.0001439045
## ... Similar to previous best
## Run 15 stress 0.0002273924
## ... Procrustes: rmse 0.000322125 max resid 0.001292199
## ... Similar to previous best
## Run 16 stress 9.993753e-05
## ... Procrustes: rmse 4.747406e-05 max resid 0.0001496008
## ... Similar to previous best
## Run 17 stress 9.9414e-05
## ... Procrustes: rmse 7.479917e-05 max resid 0.0001668549
## ... Similar to previous best
## Run 18 stress 9.94672e-05
## ... Procrustes: rmse 0.00010456 max resid 0.0002113862
## ... Similar to previous best
## Run 19 stress 9.748534e-05
## ... Procrustes: rmse 6.363321e-05 max resid 0.0001336849
## ... Similar to previous best
## Run 20 stress 9.760704e-05
## ... Procrustes: rmse 6.106859e-05 max resid 0.0001475123
## ... Similar to previous best
## *** Best solution repeated 20 times
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly)
## zero: you may have insufficient data
plot_its_trt<- plot_ordination(strat_rel_dads_its.ps,
ps_NMDS_its ,
color= "Treatment") +
#ggtitle("NMDS ordination plot - relative abundance") +
labs(color = "Treatment") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_its_trt<-plot_its_trt +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_its_trt
plot_its_inc<- plot_ordination(strat_rel_dads_its.ps,
ps_NMDS_its ,
color= "Incubation") +
#ggtitle("NMDS ordination plot - relative abundance") +
labs(color = "Incubation") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_its_inc<-plot_its_inc +
stat_ellipse(linetype =2) +
theme_bw()
plot_ellipses_its_inc
plot_its_soil<- plot_ordination(strat_rel_dads_its.ps,
ps_NMDS_its ,
color= "Soil") +
#ggtitle("NMDS ordination plot - relative abundance") +
labs(color = "Soil type") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_its_soil<-plot_its_soil +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_its_soil
Let’s proceed with ordination for the 16S dataset-NMDS
set.seed(1)
strat_rel_dads_16S.ps<-transform_sample_counts(dads_5asv_16S.ps, function(x) x / sum(x))
ps_NMDS_16S = ordinate(strat_rel_dads_16S.ps,
method="NMDS",
distance="bray")
## Run 0 stress 8.996366e-05
## Run 1 stress 8.317019e-05
## ... New best solution
## ... Procrustes: rmse 0.003178651 max resid 0.01480472
## Run 2 stress 9.962288e-05
## ... Procrustes: rmse 0.02913161 max resid 0.1366179
## Run 3 stress 9.148e-05
## ... Procrustes: rmse 0.03153671 max resid 0.1479125
## Run 4 stress 9.094561e-05
## ... Procrustes: rmse 0.02732531 max resid 0.1281317
## Run 5 stress 8.507797e-05
## ... Procrustes: rmse 0.03347368 max resid 0.1570019
## Run 6 stress 9.458946e-05
## ... Procrustes: rmse 0.01482677 max resid 0.06938939
## Run 7 stress 9.602634e-05
## ... Procrustes: rmse 0.03118334 max resid 0.1427293
## Run 8 stress 9.124274e-05
## ... Procrustes: rmse 0.008086828 max resid 0.03753532
## Run 9 stress 9.099226e-05
## ... Procrustes: rmse 0.007628467 max resid 0.03563048
## Run 10 stress 9.193224e-05
## ... Procrustes: rmse 0.03255636 max resid 0.152698
## Run 11 stress 8.232243e-05
## ... New best solution
## ... Procrustes: rmse 0.04187255 max resid 0.1900619
## Run 12 stress 9.230235e-05
## ... Procrustes: rmse 0.06475307 max resid 0.303561
## Run 13 stress 8.965179e-05
## ... Procrustes: rmse 0.0480184 max resid 0.2242121
## Run 14 stress 9.391617e-05
## ... Procrustes: rmse 0.06917108 max resid 0.3244004
## Run 15 stress 9.920624e-05
## ... Procrustes: rmse 0.0691793 max resid 0.3244382
## Run 16 stress 9.5959e-05
## ... Procrustes: rmse 0.06341095 max resid 0.2972143
## Run 17 stress 8.293653e-05
## ... Procrustes: rmse 0.06930239 max resid 0.3250205
## Run 18 stress 8.42273e-05
## ... Procrustes: rmse 0.06202637 max resid 0.2906629
## Run 19 stress 8.432815e-05
## ... Procrustes: rmse 0.06944932 max resid 0.32571
## Run 20 stress 9.080799e-05
## ... Procrustes: rmse 0.02793296 max resid 0.1292364
## *** Best solution was not repeated -- monoMDS stopping criteria:
## 20: stress < smin
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly)
## zero: you may have insufficient data
plot_16S_trt<- plot_ordination(strat_rel_dads_16S.ps,
ps_NMDS_16S ,
color= "Treatment") +
#ggtitle("NMDS ordination plot - relative abundance") +
labs(color = "Treatment") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_16S_trt<-plot_16S_trt +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_16S_trt
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
plot_16S_inc<- plot_ordination(strat_rel_dads_16S.ps,
ps_NMDS_16S ,
color= "Incubation") +
#ggtitle("NMDS ordination plot - relative abundance") +
labs(color = "Incubation") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_16S_inc<-plot_16S_inc +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_16S_inc
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
plot_16S_soil<- plot_ordination(strat_rel_dads_16S.ps,
ps_NMDS_16S ,
color= "Soil") +
#ggtitle("NMDS ordination plot - relative abundance") +
labs(color = "Soil type") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_16S_soil<-plot_16S_soil +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_16S_soil
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
Let’s proceed with ordination for the RPS10 dataset-NMDS
set.seed(1)
strat_rel_dads_rps10.ps<-transform_sample_counts(dads_5asv_rps10.ps, function(x) x / sum(x))
ps_NMDS_rps10 = ordinate(strat_rel_dads_rps10.ps,
method="NMDS",
distance="bray")
## Run 0 stress 9.646111e-05
## Run 1 stress 9.681011e-05
## ... Procrustes: rmse 0.0001489952 max resid 0.0004893922
## ... Similar to previous best
## Run 2 stress 8.725043e-05
## ... New best solution
## ... Procrustes: rmse 9.343795e-05 max resid 0.0003534458
## ... Similar to previous best
## Run 3 stress 9.218624e-05
## ... Procrustes: rmse 9.166428e-05 max resid 0.0002957858
## ... Similar to previous best
## Run 4 stress 9.131775e-05
## ... Procrustes: rmse 8.62371e-05 max resid 0.0003601615
## ... Similar to previous best
## Run 5 stress 9.687112e-05
## ... Procrustes: rmse 8.700801e-05 max resid 0.0001959933
## ... Similar to previous best
## Run 6 stress 9.390856e-05
## ... Procrustes: rmse 0.0001078327 max resid 0.0003442075
## ... Similar to previous best
## Run 7 stress 9.920391e-05
## ... Procrustes: rmse 5.200916e-05 max resid 0.000175748
## ... Similar to previous best
## Run 8 stress 9.720681e-05
## ... Procrustes: rmse 9.55066e-05 max resid 0.0002879497
## ... Similar to previous best
## Run 9 stress 8.898622e-05
## ... Procrustes: rmse 9.951896e-05 max resid 0.0002879941
## ... Similar to previous best
## Run 10 stress 8.54193e-05
## ... New best solution
## ... Procrustes: rmse 7.24128e-05 max resid 0.0002186124
## ... Similar to previous best
## Run 11 stress 9.003123e-05
## ... Procrustes: rmse 6.651538e-05 max resid 0.0002296343
## ... Similar to previous best
## Run 12 stress 8.8711e-05
## ... Procrustes: rmse 6.639742e-05 max resid 0.0002112206
## ... Similar to previous best
## Run 13 stress 9.142818e-05
## ... Procrustes: rmse 0.0001033687 max resid 0.0002963042
## ... Similar to previous best
## Run 14 stress 8.759001e-05
## ... Procrustes: rmse 0.0001045275 max resid 0.0002975806
## ... Similar to previous best
## Run 15 stress 8.993688e-05
## ... Procrustes: rmse 9.302378e-05 max resid 0.0002955304
## ... Similar to previous best
## Run 16 stress 8.753563e-05
## ... Procrustes: rmse 7.733486e-05 max resid 0.0002073648
## ... Similar to previous best
## Run 17 stress 8.508021e-05
## ... New best solution
## ... Procrustes: rmse 7.811403e-05 max resid 0.0002102992
## ... Similar to previous best
## Run 18 stress 9.858152e-05
## ... Procrustes: rmse 9.462793e-05 max resid 0.0002329584
## ... Similar to previous best
## Run 19 stress 9.620554e-05
## ... Procrustes: rmse 9.740211e-05 max resid 0.000223294
## ... Similar to previous best
## Run 20 stress 9.002458e-05
## ... Procrustes: rmse 6.780994e-05 max resid 0.0001783893
## ... Similar to previous best
## *** Best solution repeated 4 times
## Warning in metaMDS(veganifyOTU(physeq), distance, ...): stress is (nearly)
## zero: you may have insufficient data
plot_rps10_trt<- plot_ordination(strat_rel_dads_rps10.ps,
ps_NMDS_rps10 ,
color= "Treatment") +
#ggtitle("NMDS ordination plot - relative abundance") +
labs(color = "Treatment") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_rps10_trt<-plot_rps10_trt +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_rps10_trt
plot_rps10_inc<- plot_ordination(strat_rel_dads_rps10.ps,
ps_NMDS_rps10 ,
color= "Incubation") +
#ggtitle("NMDS ordination plot - relative abundance") +
labs(color = "Incubation") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_rps10_inc<-plot_rps10_inc +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_rps10_inc
plot_rps10_soil<- plot_ordination(strat_rel_dads_rps10.ps,
ps_NMDS_rps10 ,
color= "Soil") +
#ggtitle("NMDS ordination plot - relative abundance") +
labs(color = "Soil type") +
theme(axis.text=element_text(size=14), axis.title=element_text(size=14), title = element_text(size=13), legend.text = element_text(size = 18))
plot_ellipses_rps10_soil<-plot_rps10_soil +
stat_ellipse(linetype = 2) +
theme_bw()
plot_ellipses_rps10_soil
Let’s combine plots into facetted plot
create_column <- function(ord_plot) {
plot_grid(
ord_plot,
nrow = 1,
align = 'v',
axis = 'l',
rel_heights = c(1)
)
}
# Create columns for each dataset and factor
col_16S_treatment <- create_column(plot_ellipses_16S_trt)
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
col_16S_soil <- create_column(plot_ellipses_16S_soil)
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
col_16S_incubation <- create_column(plot_ellipses_16S_inc)
## Warning in MASS::cov.trob(data[, vars]): Probable convergence failure
col_its_treatment <- create_column(plot_ellipses_its_trt)
col_its_soil <- create_column(plot_ellipses_its_soil)
col_its_incubation <- create_column(plot_ellipses_its_inc)
col_rps10_treatment <- create_column(plot_ellipses_rps10_trt)
col_rps10_soil <- create_column(plot_ellipses_rps10_soil)
col_rps10_incubation <- create_column(plot_ellipses_rps10_inc)
# Create row plots with labels
row_16S <- plot_grid(col_16S_treatment, col_16S_soil, col_16S_incubation,
ncol = 3, labels = "A.", label_size = 14, hjust = -0.5, vjust = 2)
row_its <- plot_grid(col_its_treatment, col_its_soil, col_its_incubation,
ncol = 3, labels = "B.", label_size = 14, hjust = -0.5, vjust = 2)
row_rps10 <- plot_grid(col_rps10_treatment, col_rps10_soil, col_rps10_incubation,
ncol = 3, labels = "C.", label_size = 14, hjust = -0.5, vjust = 2)
# Combine all elements
facetted_figure <- plot_grid(
row_16S,
row_its,
row_rps10,
nrow = 3,
align = 'v',
axis = 'l',
rel_heights = c(1.1, 1, 1) # Slightly more height for the first row
)
facetted_figure
# Save the figure
ggsave("Figures/facetted_figure_with_ordination_NMDS.svg", facetted_figure, width = 12, height = 12)
ggsave("Figures/facetted_figure_with_ordination_NMDS.pdf", facetted_figure, width = 12, height = 12)
Software used
sessioninfo::session_info()
## ─ Session info ───────────────────────────────────────────────────────────────
## setting value
## version R version 4.4.2 (2024-10-31)
## os macOS Ventura 13.2.1
## system aarch64, darwin20
## ui X11
## language (EN)
## collate en_US.UTF-8
## ctype en_US.UTF-8
## tz America/Los_Angeles
## date 2025-02-19
## pandoc 3.2 @ /Applications/RStudio.app/Contents/Resources/app/quarto/bin/tools/aarch64/ (via rmarkdown)
## quarto 1.3.450 @ /usr/local/bin/quarto
##
## ─ Packages ───────────────────────────────────────────────────────────────────
## package * version date (UTC) lib source
## abind 1.4-8 2024-09-12 [1] CRAN (R 4.4.1)
## ade4 1.7-23 2025-02-14 [1] CRAN (R 4.4.1)
## agricolae * 1.3-7 2023-10-22 [1] CRAN (R 4.4.0)
## AlgDesign 1.2.1.1 2024-09-21 [1] CRAN (R 4.4.1)
## ape 5.8-1 2024-12-16 [1] CRAN (R 4.4.1)
## Biobase * 2.66.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## BiocGenerics * 0.52.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## BiocParallel 1.40.0 2024-11-23 [1] Bioconductor 3.20 (R 4.4.2)
## biomformat 1.34.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## Biostrings * 2.74.1 2024-12-16 [1] Bioconductor 3.20 (R 4.4.2)
## bslib 0.9.0 2025-01-30 [1] CRAN (R 4.4.1)
## cachem 1.1.0 2024-05-16 [1] CRAN (R 4.4.1)
## cli 3.6.4 2025-02-13 [1] CRAN (R 4.4.1)
## cluster 2.1.8 2024-12-11 [1] CRAN (R 4.4.1)
## codetools 0.2-20 2024-03-31 [1] CRAN (R 4.4.2)
## colorspace 2.1-1 2024-07-26 [1] CRAN (R 4.4.1)
## cowplot * 1.1.3 2024-01-22 [1] CRAN (R 4.4.0)
## crayon 1.5.3 2024-06-20 [1] CRAN (R 4.4.1)
## data.table * 1.16.4 2024-12-06 [1] CRAN (R 4.4.1)
## DelayedArray 0.32.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## DESeq2 * 1.46.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## digest 0.6.37 2024-08-19 [1] CRAN (R 4.4.1)
## doParallel * 1.0.17 2022-02-07 [1] CRAN (R 4.4.0)
## dplyr * 1.1.4 2023-11-17 [1] CRAN (R 4.4.0)
## evaluate 1.0.3 2025-01-10 [1] CRAN (R 4.4.1)
## farver 2.1.2 2024-05-13 [1] CRAN (R 4.4.1)
## fastmap 1.2.0 2024-05-15 [1] CRAN (R 4.4.1)
## foreach * 1.5.2 2022-02-02 [1] CRAN (R 4.4.0)
## furrr * 0.3.1 2022-08-15 [1] CRAN (R 4.4.0)
## future * 1.34.0 2024-07-29 [1] CRAN (R 4.4.0)
## generics 0.1.3 2022-07-05 [1] CRAN (R 4.4.1)
## GenomeInfoDb * 1.42.3 2025-01-27 [1] Bioconductor 3.20 (R 4.4.2)
## GenomeInfoDbData 1.2.13 2025-02-07 [1] Bioconductor
## GenomicRanges * 1.58.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## ggplot2 * 3.5.1 2024-04-23 [1] CRAN (R 4.4.0)
## globals 0.16.3 2024-03-08 [1] CRAN (R 4.4.1)
## glue 1.8.0 2024-09-30 [1] CRAN (R 4.4.1)
## gridExtra 2.3 2017-09-09 [1] CRAN (R 4.4.1)
## gtable 0.3.6 2024-10-25 [1] CRAN (R 4.4.1)
## hms 1.1.3 2023-03-21 [1] CRAN (R 4.4.0)
## htmltools 0.5.8.1 2024-04-04 [1] CRAN (R 4.4.1)
## httr 1.4.7 2023-08-15 [1] CRAN (R 4.4.0)
## igraph 2.1.4 2025-01-23 [1] CRAN (R 4.4.1)
## IRanges * 2.40.1 2024-12-05 [1] Bioconductor 3.20 (R 4.4.2)
## iterators * 1.0.14 2022-02-05 [1] CRAN (R 4.4.1)
## jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.4.0)
## jsonlite 1.8.9 2024-09-20 [1] CRAN (R 4.4.1)
## knitr 1.49 2024-11-08 [1] CRAN (R 4.4.1)
## labeling 0.4.3 2023-08-29 [1] CRAN (R 4.4.1)
## lattice * 0.22-6 2024-03-20 [1] CRAN (R 4.4.2)
## lifecycle 1.0.4 2023-11-07 [1] CRAN (R 4.4.1)
## listenv 0.9.1 2024-01-29 [1] CRAN (R 4.4.1)
## locfit 1.5-9.11 2025-02-03 [1] CRAN (R 4.4.1)
## magrittr * 2.0.3 2022-03-30 [1] CRAN (R 4.4.1)
## MASS 7.3-64 2025-01-04 [1] CRAN (R 4.4.1)
## Matrix 1.7-2 2025-01-23 [1] CRAN (R 4.4.1)
## MatrixGenerics * 1.18.1 2025-01-09 [1] Bioconductor 3.20 (R 4.4.2)
## matrixStats * 1.5.0 2025-01-07 [1] CRAN (R 4.4.1)
## metacoder * 0.3.8 2025-02-11 [1] CRAN (R 4.4.1)
## mgcv 1.9-1 2023-12-21 [1] CRAN (R 4.4.2)
## multtest 2.62.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## munsell 0.5.1 2024-04-01 [1] CRAN (R 4.4.1)
## nlme 3.1-167 2025-01-27 [1] CRAN (R 4.4.1)
## parallelly 1.42.0 2025-01-30 [1] CRAN (R 4.4.1)
## permute * 0.9-7 2022-01-27 [1] CRAN (R 4.4.1)
## phyloseq * 1.50.0 2024-10-29 [1] Bioconductor 3.20 (R 4.4.2)
## pillar 1.10.1 2025-01-07 [1] CRAN (R 4.4.1)
## pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.4.1)
## plyr 1.8.9 2023-10-02 [1] CRAN (R 4.4.1)
## purrr * 1.0.4 2025-02-05 [1] CRAN (R 4.4.1)
## R6 2.6.1 2025-02-15 [1] CRAN (R 4.4.1)
## ragg 1.3.3 2024-09-11 [1] CRAN (R 4.4.1)
## Rcpp 1.0.14 2025-01-12 [1] CRAN (R 4.4.1)
## readr * 2.1.5 2024-01-10 [1] CRAN (R 4.4.0)
## reshape2 1.4.4 2020-04-09 [1] CRAN (R 4.4.0)
## rhdf5 2.50.2 2025-01-09 [1] Bioconductor 3.20 (R 4.4.2)
## rhdf5filters 1.18.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## Rhdf5lib 1.28.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## rlang 1.1.5 2025-01-17 [1] CRAN (R 4.4.1)
## rmarkdown 2.29 2024-11-04 [1] CRAN (R 4.4.1)
## rstudioapi 0.17.1 2024-10-22 [1] CRAN (R 4.4.1)
## S4Arrays 1.6.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## S4Vectors * 0.44.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## sass 0.4.9 2024-03-15 [1] CRAN (R 4.4.0)
## scales 1.3.0 2023-11-28 [1] CRAN (R 4.4.0)
## sessioninfo 1.2.3 2025-02-05 [1] CRAN (R 4.4.1)
## SparseArray 1.6.1 2025-01-20 [1] Bioconductor 3.20 (R 4.4.2)
## stringi 1.8.4 2024-05-06 [1] CRAN (R 4.4.1)
## stringr * 1.5.1 2023-11-14 [1] CRAN (R 4.4.0)
## SummarizedExperiment * 1.36.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## survival 3.8-3 2024-12-17 [1] CRAN (R 4.4.1)
## svglite 2.1.3 2023-12-08 [1] CRAN (R 4.4.0)
## systemfonts 1.2.1 2025-01-20 [1] CRAN (R 4.4.1)
## textshaping 1.0.0 2025-01-20 [1] CRAN (R 4.4.1)
## tibble 3.2.1 2023-03-20 [1] CRAN (R 4.4.0)
## tidyselect 1.2.1 2024-03-11 [1] CRAN (R 4.4.0)
## tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.4.0)
## UCSC.utils 1.2.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## vctrs 0.6.5 2023-12-01 [1] CRAN (R 4.4.0)
## vegan * 2.6-10 2025-01-29 [1] CRAN (R 4.4.1)
## withr 3.0.2 2024-10-28 [1] CRAN (R 4.4.1)
## xfun 0.50 2025-01-07 [1] CRAN (R 4.4.1)
## XVector * 0.46.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
## yaml 2.3.10 2024-07-26 [1] CRAN (R 4.4.1)
## zlibbioc 1.52.0 2024-11-08 [1] Bioconductor 3.20 (R 4.4.1)
##
## [1] /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/library
## * ── Packages attached to the search path.
##
## ──────────────────────────────────────────────────────────────────────────────